%HB_KSG Replication PSRM 
%Franzese, Hays, Cook 2014
clear all;

r=100;

W=csvread('C:\Users\scott\Dropbox\Spatial Probit PSRM\PSRM Conditional Accept\Replication Folder\Application - Civil War\Data\hb_ksg_contig_W.csv');%non-rs
WNT=normw(W);


A=csvread('C:\Users\scott\Dropbox\Spatial Probit PSRM\PSRM Conditional Accept\Replication Folder\Application - Civil War\Data\hb_ksg_rep_Africa.csv');

Y=A(:,4);%Civil War Incidence 

nt = length(Y);  

XNT=A(:,[8:15]);%neighpol neighpolsq neighlgdpl lnpop polity2l polity2sq postcoldw lgdp96l 

k=size(XNT,2)+1;

TL=csvread('C:\Users\scott\Dropbox\Spatial Probit PSRM\PSRM Conditional Accept\Replication Folder\Application - Civil War\Data\hb_ksg_TL.csv');

i_nt = eye(nt); 

rand_mat1 = rand(nt,(r/2));
rand_mat2 = (1 - rand_mat1);
rand_mat = [rand_mat1, rand_mat2];


z = 1-2*Y;
Z = diag(z);


ytl = TL*Y; 

niave_ysl = WNT*Y; 


%Niave Probit 
cons = ones(nt,1); 
X_mat = [cons XNT niave_ysl]; 
[b,dev,stats] = glmfit([cons(31:nt) XNT(31:nt,:) niave_ysl(31:nt) ],Y(31:nt),'binomial','link','probit','constant','off');
beta_naive = stats.beta;
se_naive = stats.se; 
temp_naive = [beta_naive, se_naive, beta_naive./se_naive] ;
printmat(temp_naive,'Naive Results','constant neighpol neighpolsq neighlgdpl lnpop polity2l polity2sq postcoldw lgdp96l spat_lag', 'beta se t-stat')
 
p0 = beta_naive'; 

lb = [-Inf; -Inf; -Inf; -Inf; -Inf; -Inf; -Inf; -Inf; -Inf; - 1.0];
ub = [Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; Inf; 1.0]; 

options = optimset;
% Modify options setting
options = optimset(options,'Display' ,'iter');
options = optimset(options,'UseParallel','always'); 
%options = optimset(options,'FunValCheck' ,'on');
options = optimset(options,'PlotFcns' ,{  @optimplotx @optimplotfval @optimplotfirstorderopt });
options = optimset(options,'Diagnostics' ,'on');
options = optimset(options,'HessUpdate' ,'bfgs'); %dfp
%options = optimset(options,'InitialHessType' ,'scaled-identity');
options = optimset(options,'Algorithm' ,'interior-point');
options = optimset(options,'LargeScale' ,'on');
options = optimset(options,'MaxFunEvals' ,10e4);
options = optimset(options,'MaxIter' ,10e2);
options = optimset(options,'TolX', 10e-8);
[p,fval,exitflag,output,lambda,grad,hessian] = ...
fmincon(@(p)sprobit_llf_hbksg_RR_rho(p,Z,rand_mat,nt,WNT,XNT,i_nt,r),p0,[],[],[],[],lb,ub,[],options);
se = sqrt(diag(inv(hessian)));
pt = p';
vcov_hat = inv(hessian); 


temp_stpmsl = [pt, se, pt./se] ;
printmat(temp_stpmsl,'STP-MSL','constant neighpol neighpolsq neighlgdpl lnpop polity2l polity2sq postcoldw lgdp96l spat_lag', 'beta se t-stat')

% save('hb_ksg_rho_results');
